Octapartite parameterization for the spin-| square lattice Heisenberg antiferromagnet 
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We introduce a novel parameterization of the Hilbert space of a spin- 1 Heisenberg antiferromagnet 
using an octapartite description of the square lattice. This provides a systematic way to model the 
ground-state wavefunction within a truncated basis. To prove its effectiveness we study systems 
with up to 64 spins using exact diagonalization. At significantly reduced computational cost, we get 
ground-state energies within 1% of the best published values. Its non-iterative nature and lack of 
exploitation of spatial symmetries make this approach suitable for generalization to doped systems. 



Introduction: The discovery of the high-T c cuprates 
has emphasized the need to understand the physics of 
charged fermions moving in a two-dimensional (2D) back- 
ground of antiferromagnetically (AFM) coupled spins, as 
represented by a doped Hubbard modeli Since an ana- 
lytical solution seems unattainable, numerical modeling 
of large-scale, strongly correlated doped systems is of ex- 
treme importance and has been pursued in many different 
ways, each with their own advantages and disadvantages. 

For example, the powerful Monte Carlo (MC) meth- 
ods give extremely accurate information for the undoped 
AFM, providing the benchmark for the ground state (GS) 
energies and correlation functions of systems with up to 
thousands of spins,- However, they are limited by the 
sign problem if additional fermions are present. This lead 
to work on alternative optimized ways for treating the 
undoped AFM, with the latest development in MC sam- 
pling being the use of variational RVB-type ansatze.— ~— In 
the same context, density matrix renormalization group 
(DMRG) has also progressed over the years£ 

Since the interactions with the fermions arc often com- 
parable or bigger than the AFM's energy scale, most such 
AFM-specialized methods require very non-trivial mod- 
ifications to adapt to doped systems, due to technicali- 
ties of the RVB ansatz or the DMRG special boundary 
conditions. In fact, recent studies of doped systems still 
perform exact diagonalization (ED) of the full Hilbert 
space rather than trying to take advantage of either of 
these schemes^— ED studies have the huge advantage 
that they provide the GS wavefunction, from which any 
GS properties, including all correlation functions, can be 
calculated. This would make ED the numerical method 
of choice, were it not for the extreme restriction on the 
sizes of doped systems that can be currently treated^ 

A very different way to approach the doped systems 
is to start from a non-optimal AFM background descrip- 
tion, like the classical Neel state.— However, the strong 
quantum fluctuations leading to strong deviations from 
a Neel-like background are believed to be an essential 
part of the low-energy physics of the doped systems. It 
is therefore necessary to find accurate yet efficient ways 
to describe the AFM background, which can also be 
straightforwardly extended to the doped problem. 

Here we present a novel way to model the GS wavefunc- 
tion for the undoped AFM, with a high level of accuracy 
that can be systematically improved. It is so effective 



that the N = 64 wavefunction can be calculated as a nu- 
merical vector in a systematic basis with a commodity 
computer. The method is formulated in the real-space 
basis of the square lattice with periodic boundary condi- 
tions, in which other interactions are defined naturally. 
Neither translational nor point-group symmetries are ex- 
ploited. Thus, it can be used as an excellent starting 
point for studying doped models where these symmetries 
are broken, or as an efficient kernel for iterative methods, 
such as renormalization and quantum cluster theory^ 

We stress that the immediate goal is not to compete 
with the MC powerhouses in dealing with large undoped 
systems. Rather, our method provides a general insight 
in the essence of the AFM background, which can be used 
to then reduce the computational cost for larger doped 
systems with only slight loss of accuracy. Availability 
of such efficient yet accurate methods is invaluable in 
allowing extensive studies of dependence on various pa- 
rameters, to shed light on their importance and effects. 

In this paper we discuss the physical insights behind 
our method and apply it to the undoped AFM, where 
its accuracy can be easily gauged. Application to doped 
models is in progress and will be presented elsewhere^ 
With the AFM exchange taken as the energy unit, a TV- 
site AFM Heisenberg model on a square lattice reads 

H afm = ^ Si ■ Sj. (1) 

The dimension of the Hilbert space is 2^, and that of 
the commonly used S z = basis for the GS is AH/(f !) 2 . 
It is the shear size of the Hilbert space that makes this 
problem so difficult, and stimulates new approaches. 

The method: After Anderson pointed out the con- 
nection between RVB states and the projected BCS- 
type wavefunction^ the RVB framework has been im- 
plemented for an AFM torus by two approaches. The 
bottom-up route consists in the explicit optimization of 
the bond amplitudes X The top-down approach is the op- 
timization of projected wave functions^ (this showed the 
coexistence of superconducting and AFM order parame- 
ters for underdoped cuprates). So long as the AFM order 
parameter is finite, there should be another way to de- 
scribe the wavefunction without the projecting ansatz or 
the micro-managing of bond amplitudes. 

One way to specify a singlet with finite AFM order is 
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to divide the bi-partite lattice into sublattices A and B, 
and perform angular momentum addition between their 
states. If these states are classified by their total and 
z-projection quantum numbers, Sa/b G [0, ^r],m-A/B G 
[— Sa/Bi Sa/b]-, we can build 5 = singlet states: 

I ) ) = 2-/ ^g^ j 5 SA,SB 5 m A -m B \SA,m A , S B , m B ) 

The Hilbert space contains a huge number of such sin- 
glets. Since the AFM ground-state is also a singlet^ if 
we label all possible singlets of Eq. © with an index 
a, we can express the GS wave function as their linear 
combination \GS) = J2 a c a \0,0) a . We seek such an or- 
thonormal basis which is also computationally efficient. 

The staggered magnetization is a well-studied observ- 
able of the ground-state wave function. For a bi-partite 
lattice in a staggered magnetic field, this quantity is for- 
mally defined as the difference between (iha) and (mj) 
in the zero-field limit. In the absence of a staggered field, 
this quantity has alternative definitions, such as: 

m2 =N2 (E(- 1 ) M ^J = V N2 J , (3) 

where Sa/b are the total sublattice spin operators. In the 
AFM GS, m = (m 2 )^ has been extrapolated to ~ 0.3 as 
N -> oo, and increases like 1/VN to ~ 0.45 for N = 32. 16 
In terms of the total spin 5 = 5,4 + 5b, we have: 

m2 = _L(25i + 25i-5 2 ) ( 4 ) 

The ground state is a singlet, 5 = 0. It follows that the 
wave function must allocate significant weight to sectors 
with high values of Sa = Sb in order to yield such large 
values of m. In fact, for the RHS of Eq. 2] to be larger 
than the expected value, the sublattice spins Sa/b have 
to be within ^ of their maximum values. Our compu- 
tational basis takes advantage of this observation. There 
are many ways to add spins quantum mechanically, but 
not all are viable because the truncation error needs to 
be bounded systematically and the computational effort 
must be small. There are two extremes: if the original 
single-site basis of Eq. <jTJ) is used, no basis transforma- 
tion is needed but enforcing truncation based on Eq. (U) 
is costly. However, if a random basis tabulated according 
to values of Sa is used, transforming Eq. (jTJ) into the new 
basis would be costly due to the many Clebsch-Gordon 
series needed. We propose a parametrization which is a 
good compromise between these two extremes. 

Starting from the case where all ^ sublattice spins add 
to the maximum of ^ , if we want to include states with 
spin down to ^ = -j — jg, the basis must allow many 
new configurations, including ones where spins have a 
total spin of while the other spins have a total spin 
of W. This suggests that groups of ^ or fewer spins 




FIG. 1: (color online) Square lattice divided in octads (shaded 
areas). Sites positioned similarly inside octads have the same 
label. Each is surrounded by neighbors with different labels. 

must be allowed to take all possible spin quantum num- 
bers. Since a larger group would overshoot the sublat- 
tice spin below the ^ " threshold" , while a smaller group 
would introduce extra "enforcement" costs as discussed 
above, we divide the full lattice into groups of 8 sites, 
see Fig. [TJ The resulting octads repeat periodically with 
translational vectors 2a(l,±l). Each spin is identified 
by the octad it belongs to, and by its position inside the 
octad. A sublattice is composed from 4 groups of spins 
indexed with the same label, e.g. Sa = So + Si + S2 + S3. 
With this arrangement, each spin interacts with spins 
from all the 4 groups of the other sublattice, e.g. a group 
spin is always to the west /east /south/north of a spin 
in group 4/5/6/7. This partition is the minimum divi- 
sion that permits such a description and is a fundamental 
building block for the wave function in the model. Even 
though we arrived at it by looking to optimize the com- 
putation, our partition has been used in other contexts.™ 
We now straightforwardly generalize Eq. @ and write 
singlets of the total lattice as: 

7 

|0, 0, /({&})) = J2 c {s„m l} II 1$. m J (5) 

i=0 

where Si,rrii are the quantum numbers for Si, i.e. the 
total spin for group i = 0, 7. Here, C|5. ;TOi } is the product 
of the appropriate Clebsch-Gordon coefficients for the 8 
pairs of quantum numbers (Si, mi). The /({5j}) index 
on the LHS keeps track of the many ways in which these 
spins can be added to form a singlet. 

If all these singlets are kept into the computational 
basis, the GS calculation is exact. The total spin of each 
group takes values Si £ [0, an d there is no reason to 
restrict them. However, our previous discussion suggests 
that while all possible 5j values must be allowed, the 
singlets with highest GS weight are those whose total 
sublattice spin is within ^ of the maximum value. 

We therefore parameterize the singlet Hilbert space in 
terms of a completeness parameter, C s £ [0, 1], and in- 
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elude in the singlet basis only states for which Sa/b > 
^(1 — C s ). For C 8 = 1, the calculation is thus exact. For 
C s = 4, this means that the maximum number of anti- 
aligned spins in the sublattice is y . Based on the discus- 
sion for the staggered magnetization, we expect this to 
already be a good variational basis for the GS. For any 
C s =/= 1, this formulation allows all possible values of Si 
for each i = 0, 7, but restricts the ways in which they can 
combine to give the total sublattice spin. The number of 
discarded states is combinatorially large due to the large 
degeneracy of states with low sublattice spins. 

Eq. ([5]) is invariant under spin rotations, thus the anti- 
alignment is enforced here in terms of angular momentum 
addition, very different from the classical Neel picture. 
This approach also maintains the full translational sym- 
metry of the Hamiltonian. If the RVB bond-amplitude 
optimization is a bottom-up way of adding AFM order 
to a singlet^ our method is a compatible top-down ap- 
proach without the need of projecting an ansatz. 

The cost of this truncation scheme is the one-time com- 
putation of matrix elements between these basis states. 
This is acceptable because conventional diagonalization 
schemes are limited by storage, not processor speed, and 
this data is reusable for other computations that involve 
the AFM. The octapartite division of Fig.[T]is valuable for 
greatly minimizing the number of finite matrix elements 
of the Hamiltonian in the basis states of Eq. ([5]). The 
two-body interaction can be computed efficiently because 
only 2 out of 8 kets in Eq. ([5]) are changed and the delta 
functions in the Clebsch-Gordan coefficients of Eq. ([5]) al- 
low data to be aligned in a computationally efficient way. 
The overall matrix can be constructed column-by-column 
and term-by-term and is extremely sparse. The practi- 
cal implication is that the grand problem is decomposed 
into many small independent computations which can be 
massively parallelized. Table U compares the C s = h 
basis size to those of other bases commonly used by non- 
iterative methods. At C s = 4, the N — 32 system was 
solved within seconds using a generic Lanczos routine. 
Even for the 9.8 x 10 8 -size basis for N = 64, the GS vector 
can be calculated with a power method using less than 15 
matrix-vector products, each taking less than 30 minutes 
with a single cpu thread on a xeon workstation. Another 
advantage is the absence of the gapless S = 1 magnon 
excitations from the variational space, which allows the 
matrix-vector products to converge rapidly. Since these 
excitations are Goldstone modes, the required 5 = 1 ba- 
sis to describe them has roughly the same size as the 



TABLE I: Size of the C 3 = \ subspace as compared to that 
of the commonly used basis. 
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S = 


S z =0 


Full 


16 


50 


1430 


12870 


2 16 


32 


11041 


35357670 


601080390 


2 :i2 


64 


9.8 x 10 s 


5.55 x 10 16 


1.83 x 10 18 


2 64 



S = basis required for the GS, and can be generated by 
appropriately changing the Clebsch-Gordon coefficients 
in Eq. ([2]). The modeling of collective magnon excita- 
tions in the presence of carriers is model-dependent and 
will be discused elsewhere 

Results: We have performed computations for the full 
basis for TV = 16 as a benchmark for larger N values. 
For N = 32 (64) we computed up to C s = \ (C s = ±). 
While the computations with C s < 1 are variational, the 
goal here is to demonstrate that C s ~ \ is a good rule- 
of-thumb value to capture well the AFM ground-state. 

The stability of this formulation is demonstrated in 
Fig. [UJa) which shows the overlap between the GS wave 
functions for consecutive values of C s . The deviation de- 
creases exponentially with increasing C s , and for C s = j 
the overlap is « 95%. We conclude that the GS vector 
is already pointed in the correct direction even for small 
C s , and subsequent increments of C s merely result in mi- 
nor improvements. This is reinforced by the convergence 
of the GS energy demonstrated in Fig. [2Kb), which sup- 
ports the scaling law ~ e~ aCs . The error decays 
exponentially, with a rate that increases with N . The 
N = 32, 64 lines cross at C s ~ \ where the fractional er- 
ror is 10%. This exponential efficiency is amplified by 
the fact that a linear decrease in C' s leads to a combina- 
torial decrease of the basis size because of the numerous 
low-spin combinations removed from the basis. 

These GS energies are not as accurate as for established 
iterative methods, however a single-pass computation at 
C s = 1/4 already has less than 2% error. Using a lin- 
ear fit of || ~ e~ x ; , an estimate of E GS at C s = 1 is 
obtained by simple integration (see Fig. [3]). From the 
C s = 1 estimates and the N~? scaling, we extrapolate 
lmiN^^ _ —0.6671, within 0.5% of the best pub- 
lished value of -0.6692ji£ even though it is achieved at a 
significantly reduced computational cost. 

We have thus far demonstrated the exponential conver- 
gence of the GS eigenvalue and eigenvector with increas- 
ing C s , which validates our proposed variational basis. 
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FIG. 2: (color online) (a) Overlap \(GS,C S \GS,C S + -|)j 2 
between ground-state wave functions corresponding to adja- 
cent values of C s for N = 16, 32, 64. Note that for N = 64, 
the overlap is already 95% for C s = |; (b) Fractional change 
Eqs d dC S m g roun d-state energy corresponding to different 
C s cutoffs. The lines are linear fits to the data. 
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FIG. 3: (color online) GS energy computed at specific C 3 
values (symbols) and the extrapolation to C s — 1 (solid lines) 
for N — 16, 32, 64. The horizontal dashed lines are the lowest 
values from Ref. [l6l . Inset: expectation values of the non- 
commuting operators (m 2 ) and (S r ■ S r+( ^L l j). 



Because MC studies cannot produce an explicit wave- 
function, they often gauge accuracy via spin-spin corre- 
lation operators which do not commute with the Hamil- 
tonian. The insets of Fig. [3] show expectation values of 
two such operators (the staggered magnetization of Eq. [4] 
and the spin-spin correlation at maximum distance) eval- 
uated from our eigenvector for a given value of C s . Al- 
though these operators are non-commuting, their matrix 
representation is diagonal in our basis and therefore triv- 
ial to compute. The convergence of these quantities is 
essentially linear for C s < h. For greater values of C s , 
d£-(S r ■ "SV-K- -)) 1S su PP r essed faster than ~ e~ Cs so 
the post-threshold convergence is even better than for 
the GS energy. On the other hand, because m is ex- 
ploited to discard low-spin parts of the Hilbcrt space, 
its approximated value is always higher than the true 
value; convergence is only asymtotical for C s > \. Thus, 
the computational discount for the GS wavefunction at 
a small C s is achieved at the cost of lowered accuracy 



for one of many definitions of staggered magnetization, 
which does not commute with Hafm- 

Conclusions: We introduced an octa-partition of the 
square lattice, and used it to build an orthonormal sin- 
glet basis for calculating the AFM GS. We used C s to 
systematically parameterize the resulting Hilbert space, 
and proposed apriori arguments to postulate that C s ~ \ 
is the minimal value that captures the GS accurately. 
The resulting basis is combinatorially smaller than in 
any other schemes, apriorily known, and computation- 
ally very efficient. We demonstrated the stability of the 
formulation in Fig. Ufa) and identified the sources of er- 
ror. Using a single thread on a commodity computer, 
single-pass matrix computations are feasible for systems 
with up to N = 64 spins, significantly exceeding the cur- 
rent full ED record of N = 40£ 

While other established methods can yield smaller er- 
rors for the GS for larger N values, they cannot be easily 
generalized. Our formulation is non-iterative and is de- 
fined in the real space coordinate of a torus without ex- 
ploiting spatial symmetries. This allows for an economic 
generalization of this approach into iterative procedures. 

One concern about Hilbert space truncation is its ap- 
plicability to complicated models. Our formulation by it- 
self does not involve MC sampling and so it is not subject 
to the sign problem. One outstanding issue in cuprates is 
the explanation and measurement of the photoemission 
spectrum at low doping and low temperature^ Informa- 
tion about the spectrum is contained in the Krylov space 
of the model Hamiltonian acted on the AFM GSjiS thus, 
a compact description of the AFM GS in a systematic 
orthonormal space can indeed narrow the space spanned 
by calculations. Because the total spin of the AFM plus 
carriers is conserved while AFM order is expected to per- 
sist for low doping at low temperature, our method can 
be generalized to reduce the Hilbert space of the doped 
systems. The details are model-dependent and outside 
the scope of this paper; this work is progress. 
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